L <- 1
n <- 10
k <- 7
d <- L / (n - 1)
# Gerando D (Matriz que faz d2u)
solver <- matrix(0, n*n, n*n)
for (i in 1:(n * n)) {
if (i <= n || i >= (n*n - n) || i %% n <= 1) {
solver[i, i] <- 1
next
}
solver[i, i - n] <- 1
solver[i, i - 1] <- 1
solver[i, i] <- -4
solver[i, i + 1] <- 1
solver[i, i + n] <- 1
}
A <- solve(solver)Álgebra Linear Computacional - Grupo 4
Solver, o primeiro problema
Aqui, primeiramente projetamos a matriz para resolver o problema:
\[ \frac{\partial}{\partial{x}} \Big(k \frac{\partial u}{\partial x}\Big) + \frac{\partial}{\partial{y}} \Big(k \frac{\partial u}{\partial y}\Big) = f(x,y) ,\forall \text{ x, y em } \Omega : (0,1) \times (0,1) \]
\[ \\\text{condições de contorno: } u(x,y) = 0, \forall \text{ x, y em } \partial{\Omega} \]
Iniciando o Problema
Perceba que o R já possui uma função para calcular a inversa de uma matriz, a função solve.
A partir deste momento, já temos a matriz A.
Projetando a função f.
Nesse caso, escolhemos que f seja a função:
\[ f(x, y) = -\frac{e^{x}}{|y| + 1} \]
f <- matrix(0, n * n, 1)
X = seq(0, L, d)
for (i in 1:length(X)) {
for (j in 1:length(X)) {
f[i + (j - 1) * n] <- - ((exp(X[i]) / (abs(X[j]) + 1)) * ((d ** 2) / k))
}
}
solution <- A %*% fAo final, obtemos a matriz solução.
- Temos a matriz A.
- Temos o vetor força f.
Fazemos:
\[ U = Af \]
Em que \(U\) é o vetor solução.
Aplicando as Condições de Contorno
Apenas definimos a matriz cujos valores serão o do espaço (em z). Caso o ponto esteja nas bordas, o valor é zero.
x <- seq(0, L, d)
y <- seq(0, L, d)
z <- matrix(0, n, n)
for (i in 1:n) {
for (j in 1:n) {
if (i == 1 || j == 1) {
z[i, j] <- 0
} else {
z[i, j] <- solution[i + (j - 1) * n]
}
}
}Plotagem do Gráfico
Finalmente, plotamos o gráfico.
ploting_mesh <- function(solution, to_pdf=FALSE) {
fig <- plot_ly(
type = 'surface',
contours = list(
x = list(show = TRUE, start = 1.5, end = 2, size = 0.04, color = 'white'),
z = list(show = TRUE, start = 0.5, end = 0.8, size = 0.05)),
x = ~x,
y = ~y,
z = ~z)
fig <- fig %>% layout(
scene = list(
xaxis = list(nticks = 20),
zaxis = list(nticks = 4),
camera = list(eye = list(x = 0, y = -1, z = 1)),
aspectratio = list(x = .9, y = .8, z = .8)))
if (to_pdf == TRUE) {
htmlwidgets::saveWidget(widget = fig, file = "hc.html")
webshot(url = "hc.html", file = "hc.png", delay = 1, zoom = 4, vheight = 500)
} else {
return(fig)
}
}
ploting_mesh(solution)Grids e Grafos
Vamos gerar uma malha estruturada e montar a matriz de adjacência.
Código necessário
# Gerar espaço dado um intervalo
linspace <- function(from, to, nodes) {
div <- (to - from) / (nodes - 1)
return(seq(from, to, div))
}
# Gera as metrizes em x e y
generate_gen_matrices <- function(lines, columns) {
mtx <- t(matrix(linspace(1, lines, columns), columns, lines))
mty <- matrix(linspace(lines, 1, lines), lines, columns)
return(list("x"=mtx, "y"=mty))
}
# Recebe a ordem de uma matriz e devolve
# um espaço de triângulos
get_triangles <- function(lines, columns, ramdomize_vertices=FALSE) {
grid <- generate_gen_matrices(lines, columns) # Gera as matrizes
mx <- grid$x
my <- grid$y
triangles_x <- matrix(0, 1, 3) # - Eixo x das coordenadas
triangles_y <- matrix(0, 1, 3) # - Eixo y das coordenadas
triangles_v <- matrix(0, 1, 3) # - Triângulos formados
# pelos vértices
vertices <- matrix(0, lines, columns) # - Para o Grid com os
# vértices numerados
V <- rep(1:(lines*columns))
if (ramdomize_vertices == TRUE) {
set.seed(245)
V <- sample(V)
}
for (i in 1:lines) {
for (j in 1:columns) {
vertices[i, j] <- V[(i - 1) * columns + j]
}
}
# Fazendo a triangulação
for (i in (lines):(2)) {
for (j in 1:(columns-1)) {
triangles_x <- rbind(
triangles_x,
c(mx[i, j], mx[i-1, j], mx[i-1, j+1]))
triangles_x <- rbind(
triangles_x,
c(mx[i, j], mx[i, j+1], mx[i-1, j+1]))
triangles_y <- rbind(
triangles_y,
c(my[i, j], my[i-1, j], my[i-1, j+1]))
triangles_y <- rbind(
triangles_y,
c(my[i, j], my[i, j+1], my[i-1, j+1]))
triangles_v <- rbind(
triangles_v,
c(vertices[lines - i + 1, j],
vertices[lines - i + 2, j],
vertices[lines - i + 2, j+1]))
triangles_v <- rbind(
triangles_v,
c(vertices[lines - i + 1, j],
vertices[lines - i + 1, j+1],
vertices[lines - i + 2, j+1]))
}
}
# Removendo a primeira linha da cada
triangles_x <- triangles_x[2:nrow(triangles_x), ]
triangles_y <- triangles_y[2:nrow(triangles_y), ]
triangles_v <- triangles_v[2:nrow(triangles_v), ]
# Retorna toda a "situação"
return(
list("mx"=mx, "my"=my,"x_axis"=triangles_x,
"y_axis"=triangles_y, "mv"=triangles_v, "vertices"=vertices))
}
# Receive a list like:
## triangled_grid$mx == grid in x axis
## triangled_grid$my == grid in y axis
## triangled_grid$trx == coordinates in x axis
## triangled_grid$trx == coordinates in x axis
## triangled_grid$try == coordinates in y axis
## triangled_grid$mv == triangles based in order of vertices
## triangled_grid$vertices == grid with vertices
plot_mesh_grid <- function(triangled_grid) {
mx <- triangled_grid$mx
my <- triangled_grid$my
trx <- triangled_grid$x_axis
try <- triangled_grid$y_axis
mv <- triangled_grid$mv
vertices <- triangled_grid$vertices
lines <- nrow(mx)
columns <- ncol(my)
# Plot inicial, apenas o grid
plot(mx, my, xlab="x", ylab="y", main="Space Grid")
q_triangles <- length(mv) / 3
# Plotagem dos triâgulos e seus respectivos números
for (i in 1:q_triangles) {
polygon(trx[i,], try[i, ])
text(sum(trx[i,]) / 3, sum(try[i, ]) / 3, i, col="blue")
}
# Plotagem dos números dos vértices
for (i in 1:lines) {
for (j in 1:columns) {
text(mx[i, j], my[lines - i + 1, j], vertices[i, j], col="red", cex=2.1)
}
}
}
# Criação da matriz de adjacência
create_adj_matrix <- function(triangled_grid) {
lines <- nrow(triangled_grid$mx)
columns <- ncol(triangled_grid$my)
adj_matrix <- matrix(0, lines*columns, lines*columns)
q_triangles <- length(triangled$mv) / 3
for (triangle in 1:q_triangles) {
for (i in 1:3) {
for (j in i:3) {
v1 <- triangled_grid$mv[triangle, i]
v2 <- triangled_grid$mv[triangle, j]
adj_matrix[v1, v2] <- 1
adj_matrix[v2, v1] <- 1
}
}
}
return(adj_matrix)
}
# Plotagem da matriz de adjacência
plot_adj_matrix <- function(adj_matrix) {
heatmap(adj_matrix, Colv = "Rowv", Rowv = NA,
symm=TRUE, main="Matriz de Adjacência")
}Vendo funcionar
A partir daqui temos todo o necessário para gerar matrizes e fazer a triangulação.
Perceba que a ordem da matriz está em dois lugares. Isso pelo fato de termos a possibilidade de colocarmos um número diferente de linhas e colunas, para realizar mais experimentos interessantes.
# Gera uma matriz com vértices ordenados e faz a triangulação
triangled <- get_triangles(5, 5, FALSE)
plot_mesh_grid(triangled)Agora, vamos à matriz de adjacências, que foi gerada através de um gráfico de heatmap, como é visto no código principal.
Isso é pelo fato da matriz de adjacência ser composta por zeros e uns, o que facilita nesse caso.
adj <- create_adj_matrix(triangled)
plot_adj_matrix(adj)Essa é a matriz gerada.
Trocando a ordem dos vértices
Neste momento, vamos experimentar trocar a ordem dos vértices para verificar como ficará a matriz de adjacências.
triangled <- get_triangles(5, 5, TRUE)
adj <- create_adj_matrix(triangled)plot_mesh_grid(triangled)
plot_adj_matrix(adj)De fato, a matriz é gerada corretamente!
Matrizes Maiores e não Quadradas
Vamos elevar levemente o nível para tatrizes maiores:
triangled <- get_triangles(8, 8, FALSE)
adj <- create_adj_matrix(triangled)
plot_mesh_grid(triangled)
plot_adj_matrix(adj)Ok, funciona! Mas e para tamanhos de linha e coluna distintos?
Vamos testar para uma matris com 5 linhas e 6 colunas.
triangled <- get_triangles(5, 6, FALSE)
adj <- create_adj_matrix(triangled)
plot_mesh_grid(triangled)
plot_adj_matrix(adj)E agora, modificando a ordem.
triangled <- get_triangles(5, 6, TRUE)
adj <- create_adj_matrix(triangled)
plot_mesh_grid(triangled)
plot_adj_matrix(adj)Figuras mais desordenadas (Delaunay)
library("deldir")Veja a importação da biblioteca “deldir” acima, utilizada para triangulizar uma figura.
Abaixo, as coordenadas de x e y diretamente escritas no código.
x <- c(2.3,3.2,7.0,1.0,4.0,8.0)
y <- c(3.1,2.0,3.5,5.0,8.0,9.0)
dxy2 <- deldir(x,y,plot=TRUE,wl="tr")A figura acima é o plot gerado pela biblioteca “deldir”. Porém, é possível fazer uma plotagem melhor:
A partir daqui, temos alguns trabalhos:
- Plotar os vértices
- Obter o índice de cada vértice
x1 <- dxy2[["delsgs"]][["x1"]]
x2 <- dxy2[["delsgs"]][["x2"]]
y1 <- dxy2[["delsgs"]][["y1"]]
y2 <- dxy2[["delsgs"]][["y2"]]
ind1 <- dxy2[["delsgs"]][["ind1"]]
ind2 <- dxy2[["delsgs"]][["ind2"]]
len <- length(x1)
# Plotando os vértices
plot(x, y)
# Plotando as linhas
for (i in 1:len) {
polygon(c(x1[i], x2[i]), c(y1[i], y2[i]), border = "red")
}
# Matriz que terá os índices
indexes <- matrix(0, length(x), 2)
# Lendo da primeira lista
for (i in 1:len) {
indexes[ind1[i], 1] <- x1[i]
indexes[ind1[i], 2] <- y1[i]
}
# Lendo da segunda lista
for (i in 1:len) {
indexes[ind2[i], 1] <- x2[i]
indexes[ind2[i], 2] <- y2[i]
}
# Plotando os índices
for (i in 1:length(x)) {
text(indexes[i, 1], indexes[i, 2], i, col="blue", cex=1.5)
}Gerando a Matriz de Adjacências
adj <- matrix(0, length(x), length(y))
for (i in 1:len) {
adj[ind1[i], ind2[i]] <- 1
adj[ind2[i], ind1[i]] <- 1
}
heatmap(adj, Colv = "Rowv", Rowv = NA,
symm=TRUE, main="Matriz de Adjacência")Finalmente, a matriz de adjacência montada.